#include <sstream>
#include <TGraph.h>

//luminosity in fb
const int lumi=5.0;
//function to compute p-value
Double_t myfunction(Double_t *x, Double_t *par)
{
  Double_t xx=x[0]; Double_t p0=par[0]; Double_t p1=par[1];
  if(p1/p0<0.02) return TMath::Prob(2*(xx+p0),2*(p0+1));
  double bound=-3*p1;
  double boundh=3*p1;
  if(p0-3*p1<0) bound=-p0;
  TF1 func("func", "TMath::Prob(2*(x+[0]+[1]),2*([0]+1))*TMath::Gaus(x,0,[2],1)",
	   bound, boundh);
  TF1 func_norm("func_norm", "TMath::Gaus(x,0,[0],1)", bound, boundh);
  func.SetParameter(0, p0);
  func.SetParameter(1, xx);
  func.SetParameter(2, p1);
  func_norm.SetParameter(0, p1);
  return func.Integral(bound, boundh)/func_norm.Integral(bound, boundh); 
}

Double_t myfunction_inverse(Double_t *x, Double_t *par)
{
  //xx=luminosity
  //p0=signal for 1.0fb^{-1}
  //p1=background
  //p2=percentage uncertainty
  Double_t xx=x[0]; Double_t p0=par[0]; 
  Double_t p1=par[1]; Double_t p2=par[2];
  double bound=-3*p1;
  double boundh=3*p1;
  if(p0-3*p1<0) bound=-p0;

  //[0]=signal
  //[1]=QCD background
  //[2]=sys. error
  //x=integration variable

  /*
  TF1 func("func", "TMath::Prob(2*([0]+[1]+x),2*([1]+1))*TMath::Gaus(x,0,[2],1)",
	   bound, boundh);
  TF1 func_norm("func_norm", "TMath::Gaus(x,0,[0],1)", bound, boundh);
  func.SetParameter(0, xx*p0);
  func.SetParameter(1, xx*p1);
  func.SetParameter(2, xx*p1*p2);
  func_norm.SetParameter(0, xx*p1*p2);
  return func.Integral(bound, boundh)/func_norm.Integral(bound, boundh)-0.16;
  */
  return TMath::Prob(2*xx*(p0+p1), 2*(xx*p1+1))-0.16;
}


void paint_plot(TFile* fsig, TFile* fbkg, TString plot, TString png="")
{

  gStyle->SetTitleY(0.98);
  gStyle->SetTitleX(0.5);

  double beta=1.0;

  TString unique_name=fsig->GetName()+plot+"_"+png;
  TString unique_name_org=fsig->GetName()+plot+"_org";
  TString unique_name_fix=fsig->GetName()+plot+"_fix";
  
  TCanvas* c1=new TCanvas(fsig->GetName()+plot, 
			  fsig->GetName()+plot, 500, 500);
  //  c1->Divide(2,1);
  
  TH2D* psig=(TH2D*) fsig->Get(plot)->Clone();
  TH2D* pbkg=(TH2D*) fbkg->Get(plot)->Clone();

  psig->Rebin2D(4,2);
  pbkg->Rebin2D(4,2);  
  
  TH2D* gsig_fix_err=pbkg->Clone();
  
  gsig_fix_err->SetName(unique_name);
  
  TF1 calcp("calcp", myfunction, 0, 1.0, 2); 
  
  string name(fsig->GetName());
  if(name.find("Pythia")!=string::npos)
    {
      if(name.find("300")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(100);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 300 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("350")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(100);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 350 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("400")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(50);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 400 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("450")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(50);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 450 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("500")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(20);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 500 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("550")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(20);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 550 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("600")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(10);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 600 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("650")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(10);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 650 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("700")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(8);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 700 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("750")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(8);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 750 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("800")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(5);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 800 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("850")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(5);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 850 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("900")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(2);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 900 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
      else if(name.find("950")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(2);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Pythia 950 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
    }
  
  if(name.find("Herwig")!=string::npos)
    {
      if(name.find("300")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(100);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Herwig 300 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}

      else if(name.find("600")!=string::npos)
	{
	  gsig_fix_err->SetTitle("Herwig 600 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	  gsig_fix_err->SetMaximum(10);
	  gsig_fix_err->SetMinimum(0);
	}
      else if(name.find("900")!=string::npos)
	{
	  gsig_fix_err->SetMaximum(2);
	  gsig_fix_err->SetMinimum(0);
	  gsig_fix_err->SetTitle
	    ("Herwig 900 GeV RPV #tilde{g}#tilde{g} Signal Sig.");
	}
    }
  
  string plotname(plot);

  if(plotname.find("_2")!=string::npos)
    beta=2.0;

  if(plotname.find("tau31_prod")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #prod #tau_{31}^{1/4}");}  
  else if(plotname.find("tau32_prod")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #prod #tau_{32}^{1/4}");}  
  else if(plotname.find("tau21_prod")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #prod #tau_{21}^{1/4}");}  
  else if(plotname.find("tau31_sum")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1}");}  
  else if(plotname.find("tau21_sum")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #lower[-0.2]{#scale[1.1]{#Sigma}#tau_{2}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1}");}  
  else if(plotname.find("tau32_sum")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{2}");}  
  else if(plotname.find("tau31")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #tau_{31}");}
  else if(plotname.find("tau21")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #tau_{21}");}
  else if(plotname.find("tau32")!=string::npos)
    {gsig_fix_err->GetXaxis()->SetTitle("Maximum #tau_{32}");}
  
  gsig_fix_err->GetXaxis()->SetRangeUser(0,1.0);
  gsig_fix_err->GetYaxis()->SetRangeUser(0,1000); 
  gsig_fix_err->GetYaxis()->SetTitleOffset(1.5);
  gsig_fix_err->GetYaxis()->SetTitle("Minimum M_{J} cutoff");
  gsig_fix_err->GetYaxis()->SetTitleOffset(1.5);
  
  int nx=gsig_fix_err->GetNbinsX();
  int ny=gsig_fix_err->GetNbinsY();
  float increx =gsig_fix_err->GetXaxis()->GetXmax()/nx;
  float increy =gsig_fix_err->GetYaxis()->GetXmax()/ny;

  TH2D* gsig_org=gsig_fix_err->Clone();
  TH2D* gsig_fix=gsig_fix_err->Clone();
  gsig_org->SetName(unique_name_org);
  gsig_fix->SetName(unique_name_fix); 

  TH2D* gsig_org_notext=gsig_fix_err->Clone();
  TH2D* gsig_fix_notext=gsig_fix_err->Clone();
  gsig_org_notext->SetName(unique_name_org+"notext");
  gsig_fix_notext->SetName(unique_name_fix+"notext"); 


  for(int i=1; i<=nx; i++)
  for(int j=1; j<=ny; j++)
    {
      double c_nb=pbkg->GetBinContent(i, j);
      double c_nberr=pbkg->GetBinError(i, j);
      double c_ns=psig->GetBinContent(i, j);

      double c_nserr=0;
      
      c_ns=psig->IntegralAndError(0, i, j, ny+1, c_nserr);
      c_nb=pbkg->IntegralAndError(0, i, j, ny+1, c_nberr);
      
      c_ns*=lumi; c_nb*=lumi; c_nberr*=lumi; c_nserr*=lumi;
      
      double tot_err = sqrt(c_nserr*c_nserr / (c_ns*c_ns) + 
			    0.5*c_nberr*c_nberr/ (c_nb*c_nb) );
      double c_sig= c_nb>0? c_ns/sqrt(c_nb):0;
      
      
      double c_sigfunc=0;

            
      if(c_nb > 100 && c_ns/sqrt(c_nb) > 30.0)
	{
	  c_sigfunc=c_ns/sqrt(c_nb);
	}
      else
	{
	  calcp.SetParameters(c_nb, 0);
	  double c_prob=2*calcp.Eval(c_ns);
	  c_sigfunc=TMath::ErfcInverse(c_prob)*sqrt(2);
	}

      if(TMath::IsNaN(c_sigfunc) || !TMath::Finite(c_sigfunc))
	c_sigfunc=0;


      gsig_fix->SetBinContent(i,j,c_sigfunc);

      gsig_fix_notext->SetBinContent(i,j,0);
      gsig_fix_notext->SetBinError(i,j,0);

      //purely for aesthetics
      if((i)*0.1 > 0.5 && j*100 <= 100)
	{
	  gsig_fix_notext->SetBinContent
	    (i,j,gsig_fix->GetBinContent(i,j));
	}

    }
     
  gStyle->SetPaintTextFormat("3.1f");

  /*
  gsig_fix->SetMinimum(gsig_fix_err.GetMinimum());
  gsig_fix->SetMaximum(gsig_fix_err.GetMaximum());
  gsig_org->SetMinimum(gsig_fix_err.GetMinimum());
  gsig_org->SetMaximum(gsig_fix_err.GetMaximum());
  */

  /*
  c1->cd(3);
  gsig_fix_err->SetTitle("Poisson stat and error incorporated");
  gsig_fix_err->Draw("TEXT45 colz");
  */

  /*
  c1->cd(1);
  gsig_org->Draw("TEXT45 colz");
  gsig_org_notext->Draw("colz same");
  */

  c1->cd(1);
  gsig_fix->SetMarkerSize(1.6);
  gsig_fix->Draw("TEXT45 colz");
  gsig_fix_notext->Draw("col same");
  
  char lumi_text_[50];
  string filename(fbkg->GetName());
  string fname=filename.substr(9,filename.length()-5);
  
  if(filename.length()>13)
    sprintf(lumi_text_, "%3.1f\% weight cutoff #scale[0.8]{#int}L=%3.1ffb^{-1}", 
	  atof(fname.c_str()), lumi);
  //  else
  //    sprintf(lumi_text_, "#beta=1.0 20%% systematics");

  
  else
    sprintf(lumi_text_, "#beta=%3.1f       #scale[0.8]{#int}L=%3.1ffb^{-1}", beta, lumi);
  


  TPaveText *lumi_text=new TPaveText(0.5, 0, 1, 100);
  lumi_text->SetBorderSize(0);
  lumi_text->SetFillColor(0);
  lumi_text->SetFillStyle(0);
  lumi_text->AddText(lumi_text_);
  

		      /*
  TLatex* lumi_text=
    new TLatex(0.680, 77.5, lumi_text_);
  lumi_text->SetTextSize(20);
		      */
  c1->cd(1);
  lumi_text->Draw();
  c1->Update();
  c1->cd(2);
  lumi_text->Draw();
  
  c1->Update();


  c1->SaveAs(unique_name+".png");
}

void paint_dist(TFile** file, TString plot, 
		TString png="", double MJ_cut=400)
{

  gStyle->SetTitleY(0.96);
  gStyle->SetTitleX(0.4);
  TString unique_name=plot+"_"+png;
  
  TCanvas* c1=new TCanvas(plot, 
			  plot, 500, 500);
  
  TH1D *pqcd, *pp300, *pp600, *pp900, *ph300, *ph600, *ph900;

  string name(plot);

  if(MJ_cut>0)
    {
      TH2D* pqcd_=(TH2D*) file[0]->Get(plot);
      TH2D* pp300_=(TH2D*) file[1]->Get(plot);
      TH2D* pp600_=(TH2D*) file[2]->Get(plot);
      TH2D* pp900_=(TH2D*) file[3]->Get(plot);
      TH2D* ph300_=(TH2D*) file[4]->Get(plot);
      TH2D* ph600_=(TH2D*) file[5]->Get(plot);
      TH2D* ph900_=(TH2D*) file[6]->Get(plot);

      int lowbin=MJ_cut*pqcd_->GetNbinsY() /1000 +1;
      
      pqcd=pqcd_->ProjectionX(plot+"QCD", lowbin, -1);
      pqcd->SetName("crap");
      pp300=pp300_->ProjectionX(plot+"Pythia_300", lowbin, -1);
      pp600=pp600_->ProjectionX(plot+"Pythia_600", lowbin, -1);
      pp900=pp900_->ProjectionX(plot+"Pythia_900", lowbin, -1);
      ph300=ph300_->ProjectionX(plot+"Herwig_300", lowbin, -1);
      ph600=ph600_->ProjectionX(plot+"Herwig_600", lowbin, -1);
      ph900=ph900_->ProjectionX(plot+"Herwig_900", lowbin, -1);
    }

  else
    {
      pqcd=(TH1D*) file[0]->Get(plot)->Clone(plot+"QCD");
      pqcd->SetName("Crap");
      pp300=(TH1D*) file[1]->Get(plot)->Clone(plot+"Pythia_300");
      pp600=(TH1D*) file[2]->Get(plot)->Clone(plot+"Pythia_600");
      pp900=(TH1D*) file[3]->Get(plot)->Clone(plot+"Pythia_900");
      ph300=(TH1D*) file[4]->Get(plot)->Clone(plot+"Herwig_300");
      ph600=(TH1D*) file[5]->Get(plot)->Clone(plot+"Herwig_600");
      ph900=(TH1D*) file[6]->Get(plot)->Clone(plot+"Herwig_900");
    }



  pqcd->SetTitle("QCD");
  pp300->SetTitle("Pythia 300 #tilde{g}#tilde{g} #times 10");
  pp600->SetTitle("Pythia 600 #tilde{g}#tilde{g} #times 10^{2}");
  pp900->SetTitle("Pythia 900 #tilde{g}#tilde{g} #times 10^{3}");
  ph300->SetTitle("Herwig 300 #tilde{g}#tilde{g} #times 10");
  ph600->SetTitle("Herwig 600 #tilde{g}#tilde{g} #times 10^{2}");
  ph900->SetTitle("Herwig 900 #tilde{g}#tilde{g} #times 10^{3}");
  
  TH1D* c_hist[7]={pqcd, pp300, pp600, pp900, ph300, ph600, ph900};

  
  pp300->Scale(10);
  ph300->Scale(10);
  pp600->Scale(1e2);
  ph600->Scale(1e2);
  pp900->Scale(1e3);
  ph900->Scale(1e3);
  

  for(int i=0; i<7; i++)
    {
      c_hist[i]->Rebin(2);
      c_hist[i]->SetLineWidth(3);
      c_hist[i]->Scale(lumi);
      //      c_hist[i]->GetXaxis()->SetTitleSize(0.06);
      c_hist[i]->GetXaxis()->SetTitleOffset(1.1);
      if(name.find("jet_pt")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("p_{T} (GeV/c^{2})");
      else if(name.find("jet_m")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("m (GeV/c^{2})");
      else if(name.find("tau32_prod")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#prod#tau_{32}^{1/4}");
      else if(name.find("tau21_prod")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#prod#tau_{21}^{1/4}");
      else if(name.find("tau31_prod")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#prod#tau_{31}^{1/4}");
      else if(name.find("tau32_sum")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{2}");
      else if(name.find("tau21_sum")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{2}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1}");
      else if(name.find("tau31_sum")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1}");
      else if(name.find("tau21")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#tau_{21}");
      else if(name.find("tau32")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#tau_{32}");
      else if(name.find("tau31")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#tau_{31}");
    }

  

  pqcd->SetLineColor(29);  pp300->SetLineColor(kGreen+2);
  pp600->SetLineColor(6);  pp900->SetLineColor(kRed);
  ph300->SetLineColor(kGreen+2);  ph600->SetLineColor(6);
  ph900->SetLineColor(kRed);

  pqcd->SetFillColor(29);

  ph300->SetLineStyle(7);
  ph600->SetLineStyle(7);
  ph900->SetLineStyle(7);

  pqcd->Draw("lE3");
  pp300->Draw("hist l same");
  pp600->Draw("hist l same");
  pp900->Draw("hist l same");
  ph300->Draw("hist l same");
  ph600->Draw("hist l same");
  ph900->Draw("hist l same");

  c1->Update();  
  double text_xloc=0.13;
  double text_yloc=0.5;
  

  if(name.find("jet")!=string::npos ||
     name.find("tau31")!=string::npos ||
     name.find("tau21_sum")!=string::npos ||
     name.find("tau21_prod")!=string::npos ||
     name.find("tau21_2")!=string::npos )
    {
      c1->BuildLegend(0.58, 0.58, 0.87, 0.88)->SetFillStyle(0);
      text_xloc=0.58;
    }
  else  
    c1->BuildLegend(0.13, 0.58, 0.42, 0.88)->SetFillStyle(0);

  MJ_cut=400;

  char title_[200];
  if(name.find("jet_pt1")!=string::npos)
    {sprintf(title_, "Leading Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_pt2")!=string::npos)
    {sprintf(title_, "Second Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_pt3")!=string::npos)
    {sprintf(title_, "Third Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_pt4")!=string::npos)
    {sprintf(title_, "Fourth Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m1")!=string::npos)
    {sprintf(title_, "First Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m2")!=string::npos)
    {sprintf(title_, "Second Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m3")!=string::npos)
    {sprintf(title_, "Third  Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m4")!=string::npos)
    {sprintf(title_, "Fourth Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau32_prod")!=string::npos)
    {sprintf(title_, "#prod #tau_{32}^{1/4} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau31_prod")!=string::npos)
    {sprintf(title_, "#prod #tau_{31}^{1/4} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau21_prod")!=string::npos)
    {sprintf(title_, "#prod #tau_{21}^{1/4} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau32_sum")!=string::npos)
    {sprintf(title_, "#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{2} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau21_sum")!=string::npos)
    {sprintf(title_, "#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{2}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau31_sum")!=string::npos)
    {sprintf(title_, "#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau32")!=string::npos)
    {sprintf(title_, "#tau_{32} after M_{J} > %d    ", MJ_cut);}
  else if(name.find("tau31")!=string::npos)
    {sprintf(title_, "#tau_{31} after M_{J} > %d    ", MJ_cut);}
  else if(name.find("tau21")!=string::npos)
    {sprintf(title_, "#tau_{21} after M_{J} > %d    ", MJ_cut);}
  pqcd->SetTitle(title_);

  char lumi_text_[50];
  sprintf(lumi_text_, "#scale[0.8]{#int}L=%3.1ffb^{-1}", lumi);
  TLatex* lumi_text=
    new TLatex(text_xloc, text_yloc, lumi_text_);
  lumi_text->SetNDC();
  lumi_text->Draw();
  c1->Update();
  c1->SaveAs(unique_name+".png");
  /*
  pqcd;
  pp300;
  pp600;
  pp900;
  ph300;
  ph600;
  ph900;
  */

}


void paint_dist_qcd(TFile** file, TString plot, 
		TString png="", double MJ_cut=400)
{

  gStyle->SetTitleY(0.96);
  gStyle->SetTitleX(0.4);
  TString unique_name=plot+"_"+png;
  
  TCanvas* c1=new TCanvas(plot, 
			  plot, 500, 500);
  
  TH1D *pqcd, *pqcd_1, *pqcd_2, *pqcd_3, 
    *pqcd_4, *pqcd_5, *pqcd_6, *pqcd_7;

  string name(plot);

  if(MJ_cut>0)
    {
      TH2D* pqcd_=(TH2D*) file[0]->Get(plot);
      TH2D* pqcd_1_=(TH2D*) file[1]->Get(plot);
      TH2D* pqcd_2_=(TH2D*) file[2]->Get(plot);
      TH2D* pqcd_3_=(TH2D*) file[3]->Get(plot);
      TH2D* pqcd_4_=(TH2D*) file[4]->Get(plot);
      TH2D* pqcd_5_=(TH2D*) file[5]->Get(plot);
      TH2D* pqcd_6_=(TH2D*) file[6]->Get(plot);
      TH2D* pqcd_7_=(TH2D*) file[7]->Get(plot);

      int lowbin=MJ_cut*pqcd_->GetNbinsY() /1000 +1;

      if(MJ_cut<1000)
	{
	  pqcd=pqcd_->ProjectionX(plot+"QCD", lowbin, -1);
	  pqcd->SetName("crap");
	  pqcd_1=pqcd_1_->ProjectionX(plot+"1", lowbin, -1);
	  pqcd_2=pqcd_2_->ProjectionX(plot+"2", lowbin, -1);
	  pqcd_3=pqcd_3_->ProjectionX(plot+"3", lowbin, -1);
	  pqcd_4=pqcd_4_->ProjectionX(plot+"4", lowbin, -1);
	  pqcd_5=pqcd_5_->ProjectionX(plot+"5", lowbin, -1);
	  pqcd_6=pqcd_6_->ProjectionX(plot+"6", lowbin, -1);
	  pqcd_7=pqcd_7_->ProjectionX(plot+"7", lowbin, -1);
	}
      
      else
	{
	  pqcd=pqcd_->ProjectionY(plot+"QCD");
	  pqcd->SetName("crap");
	  pqcd_1=pqcd_1_->ProjectionY(plot+"1");
	  pqcd_2=pqcd_2_->ProjectionY(plot+"2");
	  pqcd_3=pqcd_3_->ProjectionY(plot+"3");
	  pqcd_4=pqcd_4_->ProjectionY(plot+"4");
	  pqcd_5=pqcd_5_->ProjectionY(plot+"5");
	  pqcd_6=pqcd_6_->ProjectionY(plot+"6");
	  pqcd_7=pqcd_7_->ProjectionY(plot+"7");
	}
    }

  else
    {
      pqcd=(TH1D*) file[0]->Get(plot)->Clone(plot+"QCD");
      pqcd->SetName("Crap");
      pqcd_1=(TH1D*) file[1]->Get(plot)->Clone(plot+"1");
      pqcd_2=(TH1D*) file[2]->Get(plot)->Clone(plot+"2");
      pqcd_3=(TH1D*) file[3]->Get(plot)->Clone(plot+"3");
      pqcd_4=(TH1D*) file[4]->Get(plot)->Clone(plot+"4");
      pqcd_5=(TH1D*) file[5]->Get(plot)->Clone(plot+"5");
      pqcd_6=(TH1D*) file[6]->Get(plot)->Clone(plot+"6");
      pqcd_7=(TH1D*) file[7]->Get(plot)->Clone(plot+"7");
    }


  pqcd->SetTitle("QCD");
  pqcd_1->SetTitle("w_{cutoff} = 8%");
  pqcd_2->SetTitle("w_{cutoff} = 5%");
  pqcd_3->SetTitle("w_{cutoff} = 2%");
  pqcd_4->SetTitle("w_{cutoff} = 1%");
  pqcd_5->SetTitle("w_{cutoff} = 0.8%");
  pqcd_6->SetTitle("w_{cutoff} = 0.5%");
  pqcd_7->SetTitle("w_{cutoff} = 0.2%");

  
  TH1D* c_hist[50]={pqcd, pqcd_1, pqcd_2, pqcd_3, pqcd_4,
		   pqcd_5, pqcd_6, pqcd_7};

  for(int i=0; i<8; i++)
    {
      // c_hist[i]->Rebin(2);
      c_hist[i]->SetLineWidth(3);
      c_hist[i]->Scale(lumi);
      c_hist[i]->SetMarkerStyle(0);
      //      c_hist[i]->GetXaxis()->SetTitleSize(0.06);
      c_hist[i]->GetXaxis()->SetTitleOffset(1.1);
      if(MJ_cut>1000)
	c_hist[i]->GetXaxis()->SetTitle("M_{J} (GeV/c^{2})");
      else if(name.find("weight")!=string::npos)
	{
	  c1->SetLogx(1);
	  c_hist[i]->GetXaxis()->SetTitle("Weight");
	}
      else if(name.find("jet_pt")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("p_{T} (GeV/c^{2})");
      else if(name.find("jet_m")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("m (GeV/c^{2})");
      else if(name.find("tau32_prod")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#prod#tau_{32}^{1/4}");
      else if(name.find("tau21_prod")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#prod#tau_{21}^{1/4}");
      else if(name.find("tau31_prod")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#prod#tau_{31}^{1/4}");
      else if(name.find("tau32_sum")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{2}");
      else if(name.find("tau21_sum")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{2}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1}");
      else if(name.find("tau31_sum")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1}");
      else if(name.find("tau21")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#tau_{21}");
      else if(name.find("tau32")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#tau_{32}");
      else if(name.find("tau31")!=string::npos)
	c_hist[i]->GetXaxis()->SetTitle("#tau_{31}");
    }

  pqcd->SetFillColor(kGray);
  pqcd->SetLineColor(kGray);  pqcd_1->SetLineColor(kBlue+2);
  pqcd_2->SetLineColor(kAzure-3); pqcd_3->SetLineColor(kTeal);
  pqcd_4->SetLineColor(kGreen);  pqcd_5->SetLineColor(kYellow);
  pqcd_6->SetLineColor(kOrange+1);  pqcd_7->SetLineColor(kRed);


  pqcd->Draw("E3 l");
  pqcd_1->Draw("hist l same");
  pqcd_2->Draw("hist l same");
  pqcd_3->Draw("hist l same");
  pqcd_4->Draw("hist l same");
  pqcd_5->Draw("hist l same");
  pqcd_6->Draw("hist l same");
  pqcd_7->Draw("hist l same");

  c1->Update();  
  double text_xloc=0.13;
  double text_yloc=0.50;

  if(MJ_cut > 1000)
    {
      c1->SetLogy(1);
    }

  if(name.find("jet")!=string::npos ||
     name.find("tau31")!=string::npos ||
     name.find("tau21_2")!=string::npos ||
     MJ_cut > 1000)
    {
      c1->BuildLegend(0.55, 0.6, 0.93, 0.88)->SetFillStyle(0);
      text_xloc=0.68;
    }
  else  
    c1->BuildLegend(0.13, 0.6, 0.48, 0.88)->SetFillStyle(0);

  if(MJ_cut<1000)
    MJ_cut=400;

  char title_[200];
  if(MJ_cut > 1000)
    {sprintf(title_, "M_{J} distributions");}
  else if(name.find("weight_bin1")!=string::npos)
    {sprintf(title_, "Weights for M_{J}#in(200,400)");}    
  else if(name.find("weight_bin2")!=string::npos)
    {sprintf(title_, "Weights for M_{J}#in(400,600)");}    
  else if(name.find("weight_bin3")!=string::npos)
    {sprintf(title_, "Weights for M_{J}#in(600,800)");}    
  else if(name.find("weight_bin4")!=string::npos)
    {sprintf(title_, "Weights for M_{J}#in(800,1000)");}    
  else if(name.find("jet_pt1")!=string::npos)
    {sprintf(title_, "Leading Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_pt2")!=string::npos)
    {sprintf(title_, "Second Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_pt3")!=string::npos)
    {sprintf(title_, "Third Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_pt4")!=string::npos)
    {sprintf(title_, "Fourth Jet p_{T} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m1")!=string::npos)
    {sprintf(title_, "First Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m2")!=string::npos)
    {sprintf(title_, "Second Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m3")!=string::npos)
    {sprintf(title_, "Third  Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("jet_m4")!=string::npos)
    {sprintf(title_, "Fourth Jet m_{j} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau32_prod")!=string::npos)
    {sprintf(title_, "#prod #tau_{32}^{1/4} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau31_prod")!=string::npos)
    {sprintf(title_, "#prod #tau_{31}^{1/4} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau21_prod")!=string::npos)
    {sprintf(title_, "#prod #tau_{21}^{1/4} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau32_sum")!=string::npos)
    {sprintf(title_, "#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{2} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau21_sum")!=string::npos)
    {sprintf(title_, "#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{2}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau31_sum")!=string::npos)
    {sprintf(title_, "#lower[-0.2]{#scale[1.1]{#Sigma}#tau_{3}}#scale[1.5]{/}#scale[1.1]{#Sigma}#tau_{1} after M_{J} > %d", MJ_cut);}
  else if(name.find("tau32")!=string::npos)
    {sprintf(title_, "#tau_{32} after M_{J} > %d    ", MJ_cut);}
  else if(name.find("tau31")!=string::npos)
    {sprintf(title_, "#tau_{31} after M_{J} > %d    ", MJ_cut);}
  else if(name.find("tau21")!=string::npos)
    {sprintf(title_, "#tau_{21} after M_{J} > %d    ", MJ_cut);}
  pqcd->SetTitle(title_);

  char lumi_text_[200];
  sprintf(lumi_text_, "#scale[0.8]{#int}L=%3.1ffb^{-1}", lumi);
  TLatex* lumi_text=
    new TLatex(text_xloc, text_yloc, lumi_text_);
  lumi_text->SetNDC(1);
  lumi_text->Draw();
  c1->Update();
  c1->SaveAs(unique_name+".png");
  /*
  pqcd;
  pqcd_1;
  pqcd_2;
  pqcd_3;
  pqcd_4;
  pqcd_5;
  */

}


TGraph* get_sig_graph(TFile** file, 
		      vector<double>& mass, int option=0, double err=0.0)
{

  //option=0 simply plot the significance
  //option=1 invert significance
  //get the luminosity for 1 sigma exclusion
  gStyle->SetTitleY(0.96);
  gStyle->SetTitleX(0.4);

  vector<double> sig;
  vector<double> sig_smear;
  vector<double> sig_err;
  vector<double> n_sig;
  vector<double> n_bkg;

  char* plotname="M_tau21_sum";
  TH2D* pbkg=(TH2D*) file[0]->Get(plotname)->Clone("QCD");
  int nx=pbkg->GetNbinsX();
  int ny=pbkg->GetNbinsY();
  TF1 calcp("calcp", myfunction, 0, 1.0, 2);
  
  pbkg->Scale(lumi);
  int taucut=0;
  if(option<0)
    taucut=nx-1;

  for(int i=0; i<mass.size(); i++)
    {
      TH2D* psig=(TH2D*) file[i+1]->Get(plotname);
      cout<<file[i+1]->GetName()<<endl;
      psig->Scale(lumi);
      double c_sig=0;
      double cn_sig=0;
      double cn_bkg=0;
      
      double c_masscut=0;
      double c_taucut=1.0;
      for(int m=taucut; m<nx; m++)
	for(int n=2; n<ny; n++)
	  {
	    double t_masscut=m*100.0;
	    double t_taucut=(n+1)*0.1;

	    double c_nbstaterr=0;
	    double c_ns=psig->Integral(0, m, n, ny+1);
	    double c_nb=pbkg->IntegralAndError(0, m, n, ny+1, c_nbstaterr);
	    if(option<0)
	      {
		c_ns=psig->Integral(0, nx+1, n, ny+1);
		c_nb=pbkg->Integral(0, nx+1, n, ny+1);
	      }

	    if(c_nbstaterr / c_nb > 0.5)
	      continue;
	    double c_sigfunc=0;
	    
	    double c_nberr=err;

	    
	    if(c_nb > 100 && c_ns/sqrt(c_nb) > 30.0)
	      {
		c_sigfunc=c_ns/sqrt(c_nb);
	      }
	    else		    
	      {
		calcp.SetParameters(c_nb, 0);
		double c_prob=2*calcp.Eval(c_ns);
		c_sigfunc=TMath::ErfcInverse(c_prob)*sqrt(2);
	      }


	    /*	    
	    calcp.SetParameters(c_nb, c_nberr*c_nb);
	    double c_prob=2*calcp.Eval(c_ns);
	    
	    if(c_prob<1e-150 || TMath::IsNaN(c_prob))
	      {
		calcp.SetParameters(c_nb/10, c_nberr*c_nb/10);
		c_prob=2*calcp.Eval(c_ns/10);
		if(TMath::IsNaN(c_prob) || c_prob < 1e-100)
		  {
		    calcp.SetParameters(c_nb/500, c_nberr*c_nb/500);
		    c_prob=2*calcp.Eval(c_ns/500);
		    if(TMath::IsNaN(c_prob) || c_prob < 1e-100)
		      c_sigfunc=0;
		    else
		      {
			//cout<<"Prob after 500: "<<c_prob<<endl;
			c_sigfunc=sqrt(500)*TMath::ErfcInverse(c_prob)*sqrt(2);
		      }
		  }
		else
		  {
		    //cout<<"Prob after 10: "<<c_prob<<endl;
		    c_sigfunc=sqrt(10)*TMath::ErfcInverse(c_prob)*sqrt(2);
		  }
	      }
	    else
	      {
		//cout<<"Prob: "<<c_prob<<endl;
		c_sigfunc=TMath::ErfcInverse(c_prob)*sqrt(2);
	      }
	    */

	    if(c_sigfunc > c_sig)
	      {
		c_sig=c_sigfunc;
		c_masscut=t_masscut;
		c_taucut=t_taucut;
		cn_sig=c_ns;
		cn_bkg=c_nb;
	      }
	  }
      sig.push_back(c_sig);
      n_sig.push_back(cn_sig);
      n_bkg.push_back(cn_bkg);
    }
  

  if(option!=0)
    {
      TF1 calcinv("calcinv", myfunction_inverse, 0, 50.0, 3);
      for(int i=0; i<sig.size(); i++)
	{
	  calcinv.SetParameters(n_sig[i]/lumi, n_bkg[i]/lumi, 0);
	  ROOT::Math::WrappedTF1 wf1(calcinv);	  
	  // Create the Integrator
	  ROOT::Math::BrentRootFinder brf;
	  
	  brf.SetFunction( wf1, 0, 20);
	  brf.Solve();

	  sig[i]=brf.Root();
	}
    }

  TGraph* plot=new TGraph(mass.size(), &mass[0], &sig[0]);
  plot->SetName("sig_mass");
  plot->SetTitle("Signal Significance vs. Mass");
  
  return plot;


}


void bin(bool redo=false)
{
  //gStyle->SetCanvasColor(-1);
  //gStyle->SetPadColor(-1);
gStyle->SetFrameFillColor(-1);
//gStyle->SetHistFillColor(-1);
//gStyle->SetTitleFillColor(-1);
//gStyle->SetFillColor(-1);
//gStyle->SetFillStyle(4000);
gStyle->SetStatStyle(0);
//gStyle->SetTitleStyle(0);
gStyle->SetCanvasBorderSize(0);
gStyle->SetFrameBorderSize(0);
gStyle->SetLegendBorderSize(0);
gStyle->SetStatBorderSize(0);
gStyle->SetTitleBorderSize(0); 

  int NRGBs = 8, NCont = 99;
  gStyle->SetNumberContours(99);
  Double_t stops[8] = { 0.00, 0.10, 0.25, 0.45, 0.55, 0.75, 0.90, 1.00};
  Double_t red[8]   = { 1.00, 0.00, 0.00, 0.00, 0.97, 0.97, 1.0, 1.0};
  Double_t green[8] = { 1.00, 0.97, 0.30, 0.40, 0.97, 0.5, 0.3, 0.00};
  Double_t blue[8]  = { 1.00, 0.97, 0.97, 0.00, 0.00, 0.00, 0.0, 0.00};
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);

  gStyle->SetOptStat(0);
  TFile* qcd=new TFile("QCD.root");
  TFile* p_300=new TFile("Pythia_glu300_ntuple.root");
  TFile* p_600=new TFile("Pythia_glu600_ntuple.root");
  TFile* p_900=new TFile("Pythia_glu900_ntuple.root");
  TFile* h_300=new TFile("Herwig_glu300_ntuple.root");
  TFile* h_600=new TFile("Herwig_glu600_ntuple.root");
  TFile* h_900=new TFile("Herwig_glu900_ntuple.root");

  TFile* c_300=new TFile("gluino_300_ntuple.root");
  TFile* c_400=new TFile("gluino_400_ntuple.root");
  TFile* c_600=new TFile("gluino_600_ntuple.root");
  TFile* c_650=new TFile("gluino_650_ntuple.root");
  TFile* c_850=new TFile("gluino_850_ntuple.root");


  if(redo){
 
    
    /*
    ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_1.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_2.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_5.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_8.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.1.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.2.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.5.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.8.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.01.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.02.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.05.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.08.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.001.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.002.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.005.root");
  ((TTree*)qcd->Get("Ntuple"))->Process("Ntuple.C++", "plot_QCD_0.008.root");
  
  ((TTree*)p_300->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_300.root");
  ((TTree*)p_600->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_600.root");
  ((TTree*)p_900->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_900.root");
  ((TTree*)h_300->Get("Ntuple"))->Process("Ntuple.C++", "plot_Herwig_300.root");
  ((TTree*)h_600->Get("Ntuple"))->Process("Ntuple.C++", "plot_Herwig_600.root");
  ((TTree*)h_900->Get("Ntuple"))->Process("Ntuple.C++", "plot_Herwig_900.root");
    */
  ((TTree*)c_300->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_300_new.root");
  ((TTree*)c_400->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_400_new.root");
  //  ((TTree*)c_600->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_600_new.root");
  ((TTree*)c_650->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_650_new.root");
  //  ((TTree*)c_850->Get("Ntuple"))->Process("Ntuple.C++", "plot_Pythia_850_new.root");
  }

  TFile* fc_300=new TFile("plot_Pythia_300_new.root");
  TFile* fc_400=new TFile("plot_Pythia_400_new.root");
  TFile* fc_650=new TFile("plot_Pythia_650_new.root");
  //  TFile* c_300=new TFile("plot_Pythia_300_new.root");
  //  TFile* c_300=new TFile("plot_Pythia_300_new.root");

  
  TFile* fqcd=new TFile("plot_QCD.root");
  TFile* fp_300=new TFile("plot_Pythia_300.root");
  TFile* fp_600=new TFile("plot_Pythia_600.root");
  TFile* fp_900=new TFile("plot_Pythia_900.root");
  TFile* fh_300=new TFile("plot_Herwig_300.root");
  TFile* fh_600=new TFile("plot_Herwig_600.root");
  TFile* fh_900=new TFile("plot_Herwig_900.root");

  

  TFile* filelist[7]={fqcd, fp_300, fp_600, fp_900, 
		      fh_300, fh_600, fh_900};
  
  TFile* fqcd_1=new TFile("plot_QCD_8.root");
  TFile* fqcd_2=new TFile("plot_QCD_5.root");
  TFile* fqcd_3=new TFile("plot_QCD_2.root");
  TFile* fqcd_4=new TFile("plot_QCD_1.root");
  TFile* fqcd_5=new TFile("plot_QCD_0.8.root");
  TFile* fqcd_6=new TFile("plot_QCD_0.5.root");
  TFile* fqcd_7=new TFile("plot_QCD_0.2.root");

  TFile* filelist_qcd[20]={fqcd, fqcd_1, fqcd_2,
			  fqcd_3, fqcd_4, fqcd_5, fqcd_6, fqcd_7};
  

  //  paint_dist_qcd(filelist_qcd, "M_tau32_2", "tau32_weight_1", 400);
  //  paint_dist_qcd(filelist_qcd, "M_tau21_1", "Mass_weight_1", 2000);

  /*
  paint_dist_qcd(filelist_qcd, "M_tau21_1", "Mass_weight_1", 2000);
  
  paint_dist_qcd(filelist_qcd, "M_tau21_2", "tau21_weight_1", 400);
  paint_dist_qcd(filelist_qcd, "M_tau32_2", "tau32_weight_1", 400);
  paint_dist_qcd(filelist_qcd, "M_tau31_2", "tau31_weight_1", 400);
  */

  /*
  paint_dist_qcd(filelist_qcd, "weight_bin1", "_weight", -1);
  paint_dist_qcd(filelist_qcd, "weight_bin2", "_weight", -1);
  paint_dist_qcd(filelist_qcd, "weight_bin3", "_weight", -1);
  paint_dist_qcd(filelist_qcd, "weight_bin4", "_weight", -1);
  */

  // paint_plot(fp_600, fqcd, "M_tau32_2", "_poisson_0.008");
  
  
  paint_plot(fc_300, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fc_300, fqcd, "M_tau21_sum", "8-9_2D");

  /*
  paint_plot(fc_400, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fc_400, fqcd, "M_tau21_sum", "8-9_2D");

  paint_plot(fp_600, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fp_600, fqcd, "M_tau21_sum", "8-9_2D");

  paint_plot(fc_650, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fc_650, fqcd, "M_tau21_sum", "8-9_2D");

  paint_plot(fp_900, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fp_900, fqcd, "M_tau21_sum", "8-9_2D");
  
  paint_plot(fp_300, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fp_300, fqcd, "M_tau21_sum", "8-9_2D");

  paint_plot(fp_300, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fp_300, fqcd, "M_tau21_sum", "8-9_2D");
  */

  /*  
  paint_plot(fp_600, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fp_600, fqcd, "M_tau21_sum", "8-9_2D");

  paint_plot(fp_900, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fp_900, fqcd, "M_tau21_sum", "8-9_2D");


  paint_plot(fh_300, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fh_300, fqcd, "M_tau21_sum", "8-9_2D");

  paint_plot(fh_600, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fh_600, fqcd, "M_tau21_sum", "8-9_2D");
  
  paint_plot(fh_900, fqcd, "M_tau32_2", "8-9_2D");
  paint_plot(fh_900, fqcd, "M_tau21_sum", "8-9_2D");
  */

  /*  
  paint_plot(fp_300, fqcd, "M_tau31_sum", "8-9_2D");
  paint_plot(fp_600, fqcd, "M_tau31_sum", "8-9_2D");
  paint_plot(fp_900, fqcd, "M_tau31_sum", "8-9_2D");
  paint_plot(fh_300, fqcd, "M_tau31_sum", "8-9_2D");
  paint_plot(fh_600, fqcd, "M_tau31_sum", "8-9_2D");
  paint_plot(fh_900, fqcd, "M_tau31_sum", "8-9_2D");

  paint_plot(fp_300, fqcd, "M_tau32_sum", "8-9_2D");
  paint_plot(fp_600, fqcd, "M_tau32_sum", "8-9_2D");
  paint_plot(fp_900, fqcd, "M_tau32_sum", "8-9_2D");
  paint_plot(fh_300, fqcd, "M_tau32_sum", "8-9_2D");
  paint_plot(fh_600, fqcd, "M_tau32_sum", "8-9_2D");
  paint_plot(fh_900, fqcd, "M_tau32_sum", "8-9_2D");
  */

  //  paint_plot(fp_600, fqcd, "M_tau32_2", "_poisson");  
  /*
  paint_plot(fp_600, fqcd_1, "M_tau32_2", "_poisson1");
  paint_plot(fp_600, fqcd_2, "M_tau32_2", "_poisson2");
  paint_plot(fp_600, fqcd_3, "M_tau32_2", "_poisson3");
  paint_plot(fp_600, fqcd_4, "M_tau32_2", "_poisson4");
  paint_plot(fp_600, fqcd_5, "M_tau32_2", "_poisson5");
  paint_plot(fp_600, fqcd_6, "M_tau32_2", "_poisson6");
  paint_plot(fp_600, fqcd_7, "M_tau32_2", "_poisson7");
  */

  /*
  paint_plot(fp_600, fqcd, "M_tau21_sum", "_poisson");
  paint_plot(fp_600, fqcd, "M_tau31_sum", "_poisson");

  paint_plot(fh_600, fqcd, "M_tau21_sum", "_poisson");
  paint_plot(fh_600, fqcd, "M_tau31_sum", "_poisson");
  */


  /*  
  paint_dist(filelist, "M_tau32_1", "8-7");
  paint_dist(filelist, "M_tau31_1", "8-7");
  paint_dist(filelist, "M_tau21_1", "8-7");

  paint_dist(filelist, "M_tau21_sum", "8-7");
  paint_dist(filelist, "M_tau32_sum", "8-7");
  paint_dist(filelist, "M_tau31_sum", "8-7");

  paint_dist(filelist, "M_tau21_prod", "8-7");
  paint_dist(filelist, "M_tau32_prod", "8-7");
  paint_dist(filelist, "M_tau31_prod", "8-7");
  */

  TFile** mass_file_list[49]={fqcd, fc_300, fc_400, fp_600,
			     fc_650, fp_900};

  vector<double> mass_list;
  mass_list.push_back(300);
  mass_list.push_back(400);
  mass_list.push_back(600);
  mass_list.push_back(650);
  //  mass_list.push_back(850);
  mass_list.push_back(900);
  
  
  TGraph* sig_plot=get_sig_graph(mass_file_list, mass_list, -1);
  TGraph* sig_plot1=get_sig_graph(mass_file_list, mass_list, 1);
  //TGraph* sig_plot1=get_sig_graph(mass_file_list, mass_list, 0, 0.1);
  //  TGraph* sig_plot2=get_sig_graph(mass_file_list, mass_list, 0, 0.2);
  //  TGraph* sig_plot3=get_sig_graph(mass_file_list, mass_list, 0, 0.3);
  

  TCanvas* a=new TCanvas("a", "a", 500, 500);
  sig_plot->Draw("AL*");
  sig_plot1->Draw("L* same");
  //  sig_plot1->Draw("same L*");
  //  sig_plot1->Draw("AL* same");
  //  sig_plot2->Draw("AL* same");
  //  sig_plot3->Draw("AL* same");
  


  /*
  paint_dist(filelist, "M_tau21_1", "_MJ400");
  paint_dist(filelist, "M_tau31_1", "_MJ400");  
  paint_dist(filelist, "M_tau32_2", "_MJ400");
  paint_dist(filelist, "M_tau21_2", "_MJ400");
  paint_dist(filelist, "M_tau31_2", "_MJ400");
  paint_dist(filelist, "M_tau31_prod", "_MJ400");
  paint_dist(filelist, "M_tau32_prod", "_MJ400");
  paint_dist(filelist, "M_tau21_prod", "_MJ400");
  paint_dist(filelist, "M_tau31_sum", "_MJ400");
  paint_dist(filelist, "M_tau32_sum", "_MJ400");
  paint_dist(filelist, "M_tau21_sum", "_MJ400");
  

  
  paint_dist(filelist, "jet_pt1", "_pt1", -1);
  paint_dist(filelist, "jet_pt2", "_pt2", -1);
  paint_dist(filelist, "jet_pt3", "_pt3", -1);
  paint_dist(filelist, "jet_pt4", "_pt4", -1);
  paint_dist(filelist, "jet_m1", "_m1", -1);
  paint_dist(filelist, "jet_m2", "_m2", -1);
  paint_dist(filelist, "jet_m3", "_m3", -1);
  paint_dist(filelist, "jet_m4", "_m4", -1);
  */
  

  
}


